The "black holes" below are simulated using different planar dynamical systems. These obviously do not take into account things like
time dilation, 3D coordinates, or the curvature of space-time. However, they are much faster to simulate, and capture the attracting and (maybe fictitious) spiral nature of a black hole.
We can change the system's parameters to get differently behaving "black holes". But how do we go about determining which ODEs to use? If you are interested in how I derived these equations to exhibit specific behavior, check under the simulation for an explanation. For now, consider the generalized planar system:
The goal is to obviously have a stable equilibrium, or possibly even two. I want to have several different types of "black holes". For example, one that attracts things in a somewhat straight line to the origin,
one with a spiraling to the center, and one with a psuedo event horizon. The first should be easily done with a system whose linearized version has negative real eigenvalues.
The second one can be done if the linearized version has complex conjugate eigenvalues with negative real parts. The final one will be the most difficult, we will use a limit cycle around
an equilibrium to generate the sort of event horizon effect. We do not have to worry about equilibria that are not hyperbolic right now, as we can force the eigenvalue to be non-zero. To get started, the Jacobian matrix is:
With characteristic polynomial:
\[\lambda^2 - \tau \lambda + \Delta = 0\]
where \(\tau = \frac{\partial f}{\partial x} + \frac{\partial g}{\partial y}\) and \(\Delta = \frac{\partial f}{\partial x} \cdot \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \cdot \frac{\partial g}{\partial x}\)
If we just want to consider real roots, then a sufficient condition for an attracting equilibrium is that
\[\tau < 0 \text{ and } \Delta > 0\]
This corresponds to
\[\frac{\partial f}{\partial x} < -\frac{\partial g}{\partial y} \text{ and } \frac{\partial f}{\partial x} \frac{\partial g}{\partial y} > \frac{\partial f}{\partial y} \frac{\partial g}{\partial x}\]
At the equilibrium. For easy initial analysis, let's try to make the equilibrium the origin. Then any mixed polynomial with no constant term will have a root at the origin.
\[f(x,y) = \alpha x^2 + \beta y^2 + \gamma xy + \omega x + \eta y\]
\[g(x,y) = \delta x^2 + \chi y^2 + \xi xy + \theta x + \sigma y\]
Where all the greek letters are real valued parameters. Then,
\[\frac{\partial f}{\partial x} = 2\alpha x + \gamma y + \omega\]
\[\frac{\partial f}{\partial y} = 2\beta y + \gamma x + \eta\]
\[\frac{\partial g}{\partial x} = 2\delta x + \xi y + \theta\]
\[\frac{\partial g}{\partial y} = 2\chi y + \xi x + \sigma\]
Evaluated at \((0,0)\) gives:
\[\frac{\partial f}{\partial x} = \omega\]
\[\frac{\partial f}{\partial y} = \eta\]
\[\frac{\partial g}{\partial x} = \theta\]
\[\frac{\partial g}{\partial y} = \sigma\]
So we need:
\[\omega < -\theta\]
and
\[\omega \sigma > \eta \theta\]
Using this criteria, we can generate different systems with slightly different behaviors. For simplicity, let \(\theta = 0\). This causes \(\omega\) and \(\sigma\) to be negative, since we need \(\omega \sigma > 0\).